Predation strongly limits demography of a keystone migratory herbivore in a recovering transfrontier ecosystem

Abstract Large herbivore migrations are imperiled globally; however the factors limiting a population across its migratory range are typically poorly understood. Zambia's Greater Liuwa Ecosystem (GLE) contains one of the largest remaining blue wildebeest (Connochaetes taurinus taurinus) migrations, yet the population structure, vital rates, and limiting factors are virtually unknown. We conducted a long‐term demographic study of GLE wildebeest from 2012 to 2019 of 107 collared adult females and their calves, 7352 herd observations, 12 aerial population surveys, and concurrent carnivore studies. We applied methods of vital rate estimation and survival analysis within a Bayesian estimation framework. From herd composition observations, we estimated rates of fecundity, first‐year survival, and recruitment as 68%, 56%, and 38% respectively, with pronounced interannual variation. Similar rates were estimated from calf‐detections with collared cows. Adult survival rates declined steadily from 91% at age 2 years to 61% at age 10 years thereafter dropping more sharply to 2% at age 16 years. Predation, particularly by spotted hyena, was the predominant cause of death for all wildebeest ages and focused on older animals. Starvation only accounted for 0.8% of all unbiased known natural causes of death. Mortality risk differed substantially between wet and dry season ranges, reflecting strong spatio‐temporal differences in habitat and predator densities. There was substantial evidence that mortality risk to adults was 27% higher in the wet season, and strong evidence that it was 45% higher in the migratory range where predator density was highest. The estimated vital rates were internally consistent, predicting a stable population trajectory consistent with aerial estimates. From essentially zero knowledge of GLE wildebeest dynamics, this work provides vital rates, age structure, limiting factors, and a plausible mechanism for the migratory tendency, and a robust model‐based foundation to evaluate the effects of potential restrictions in migratory range, climate change, predator–prey dynamics, and poaching.


| INTRODUC TI ON
Large ungulates play a critical role in the functioning of ecosystems through their impacts on vegetation structure, diversity, and nutrient cycling, and as prey for top carnivores, yet they are some of the most imperiled species worldwide (Ripple et al., 2015). Ungulates are most abundant in migratory populations, and seasonal migration is thought to mitigate the regulatory effects of food limitation and predation Hopcraft et al., 2015). While previous research has focused on spatio-temporal variation in food resources and predation, the connection between migratory behavior and population limitation is unknown for most large herbivore migrations (Bolger et al., 2008). Rapid human-induced ecological change has made the long-term prospects for large ungulates in general-and migratory populations in particular-very poor given that they are wide-ranging, often predictable in their migratory routes and timing, and rarely have their entire migratory range protected (Berger, 2004;Hopcraft et al., 2015;Ripple et al., 2015).
Such characteristics make these populations extremely sensitive to human impacts in the form of overhunting, poaching, human barriers such as fencing, and habitat loss from agricultural conversion and other practices (Bolger et al., 2008), particularly in the face of climate change (Durant et al., 2015). In Africa the current decline of migratory herbivores is most pronounced (Ripple et al., 2015) and exemplified by the iconic migratory ungulate, the blue wildebeest (Connochaetes taurinus). Once widespread throughout savannah Africa the wildebeest has experienced range-wide declines, specifically when migratory routes have been stopped or altered, resulting in rapid population collapse (Berger, 2004;Bolger et al., 2008;Morrison et al., 2016). While these declines have been well documented, the factors driving demography of migratory wildebeest populations-and their subsequent declines-are not well understood.
Long-lived mammals such as ungulates exhibit age and sexdependent vital rates (Eberhardt, 2002) and these rates can be impacted by an array of ecological factors. These broadly include resource limitation from forage and minerals as well as predation and predation risk effects (e.g., changes in behavior, habitat selection, nutrition and reproduction as a result of a predator being on the landscape, Creel et al., 2019), all of which can interact and in turn can be affected by and interact with human activities. To better conserve rapidly declining migratory populations, Bolger et al. (2008) emphasized an integrated approach in part aimed at understanding population limitation in all phases of the migratory cycle and the demographic consequences of human disruption to these migrations. Despite blue wildebeest being very well-studied through the seminal work on the Serengeti-Mara populations (Hopcraft et al., 2015;Mduma et al., 1999;Sinclair & Arcese, 1995a, 1995bSinclair & Norton-Griffiths, 1979) few intensive, long-term, and individual-based demographic studies exist on this keystone species to evaluate population limitation in this context (Bolger et al., 2008). And more such long-term demographic studies are needed on large ungulates in general (Gaillard, 2013;Schradin & Hayes, 2017).
Western Zambia's Liuwa Plain National Park (LPNP) and the surrounding area comprises the Greater Liuwa Ecosystem (GLE), which currently houses what is thought to be the second-largest remaining wildebeest migration in Africa, and the largest population of common blue wildebeest (Estes, 2014). It is part of a larger ecosystem spanning much of northeastern Angola and comprising the Liuwa-Mussuma Transfrontier Conservation Area (LMTFCA, M'soka et al., 2016). Historically wildebeest are thought to have seasonally migrated long distances between both countries prior to being decimated in the decades-long Angolan civil war , and the boundaries of the LMTFCA were in part designed to protect the historical migration route and patterns, although this migration is largely undocumented (Harris et al., 2009) However, the ecological and anthropogenic factors limiting this migratory population across its range are not well understood, and this has significant implications for conservation and management of both this keystone species and ecosystem. Wildebeest are experiencing an array of limiting demographic impacts in the form of predation from a recovering carnivore population  2 | ME THODS

| Study area
The LMTFCA is located on the boundary of Western Zambia and Angola, and our study area comprised the Greater Liuwa Declared a National Park in 1972, LPNP has been a conservation area since the 1880s, when it was declared by the Lozi King as a royal hunting ground. Following park designation, settlements were still permitted and the human population within and around LPNP numbered approximately 17,000 people during this study (African Parks, 2018). Largely as a result of a decades-long civil war in neighboring Angola (1975Angola ( -2002 the wildlife populations in the GLE were decimated, with many species extirpated or reduced to severely low numbers (M'soka et al., 2016). In 2003, African Parks signed a 20year public-private partnership with the DNPW with the goal of restoring the GLE ecosystem and making it a profitable wildlife-based economy. The GLE wildebeest population is thought to have been transboundary across the border with Angola. Currently it exists only within Zambia, migrating seasonally between a rainy season range generally in the southern portion of LPNP and a dry season range in northern LPNP and the adjacent UWZGMA (M'soka et al., 2017). Liuwa also contains an unknown but large number of oribi (Ourebia ourebi) and approximately 5000 migratory zebra (Equus burchelli), as well as smaller numbers of red lechwe (Kobus leche leche), tsessebe (Damaliscus lunatus lunatus), reedbuck (Redunca arundium), common duiker (Sylvicapra grimmia), buffalo (Syncerus caffer), eland (Tragelaphus oryx), roan (Hippotragus equinus), and sitatunga (Tragelaphus spekii) (Viljoen, 2009(Viljoen, , 2011(Viljoen, , 2013. A 1970 aerial census counted 12,691 wildebeest and estimated a population as high as 16,000 (Estes & East, 2009). More recently, 10 aerial surveys between 2001 and 2018 estimated wildebeest and other herbivores, with best estimates ranging between 23,000 and 46,000 and wide confidence intervals where these were quantified (African Parks, 2018;Viljoen, 2013). Almost the entire population is located in the southeast portion of the range from December to April; migration to the northwest occurs in May and June and is complete by July. Migration to the southeast occurs in October and November. Calving occurs around the first of October, and evidence of weaning (i.e., spatial dissociation of calves from their mothers) begins around July.

F I G U R E 1 Study area in and around
Liuwa plain National Park, Northwest Zambia. Concentric circles indicate distance from a point beyond the southeastern extreme of wildebeest range, as a basis for a one-dimensional measure of location along the migratory corridor.

| Capture and collaring
We focused our study around adult female wildebeest and utilized radio-collars for the basis of this work. We initiated collaring in May 2012 when a study population of 45 adult cows were fitted with Telonics VHF collars with a 5-year battery life. Because losses to predation and poaching occurred throughout the study, we attempted to maintain a study population of 30-50 adult cows year-round and thus conducted additional collar deployments of 6-15 animals from October to December in ensuing years, following calving when most of the population was congregated in the south. A total of 117 collars were deployed on 107 cows, including 95 VHF collars, 18 GPS collars, four satellite/GPS collars, and 10 recollarings of previously collared animals. Animals were immobilized with cartridge-fired dart rifles projected from vehicle or helicopter in the initial capture, and then by vehicle thereafter. Immobilized wildebeest were monitored for temperature and respiration and fitted with collars while blood, tissue, hair, and fecal samples were obtained from each animal.
Measurements of shoulder height and body length were performed, as well as measurements of tooth wear on the I1 and I2 incisor teeth to estimate age (Christianson et al., 2018). When possible, the animal was physically checked for pregnancy. The immobilization was then reversed, and the animal was monitored for recovery. All procedures followed guidelines established by the Department of Veterinary Services (Zambia) and the Department of National Parks and Wildlife (Zambia) and a Zambian-registered veterinarian conducted and oversaw all operations.

| Aging
Given that large ungulates exhibit age-specific rates of survival and reproduction (Eberhardt, 2002), it was important to have an estimate of age for each cow. Thus, we developed a novel aging method based on arcade breadth, labio-lingual width, and crown-height measurements of the first incisors that correlated well with cementum annuli analyses from wildebeest that had been predated during the course of the carnivore work and their incisors removed for analyses (Christianson et al., 2018). Consequently, when wildebeest were immobilized a trained team collected these measurements on both first incisors to provide the basis for an estimate of the cow's age using an algorithm detailed in the Appendix A.

| Field methods for demographic data collection
Beginning in September 2012 we attempted to visually locate each collared cow on the ground at least once every 2 weeks. Field teams covered the study area daily by motorbike or vehicle to radio-track collared cows, and periodic aerial tracking flights were also utilized when possible using an ultralight aircraft (Foxbat, Aeroprakt, Ukraine). Given that the population was migratory and unstudied, it was difficult to implement a stratified random sampling procedure of collared cows. Thus, we employed an intensive search effort throughout the study area to locate every animal regularly, logging approximately 800-1100 person days in the field each year.
We not only collected survival and reproduction data on these collared cows but also used these locations to collect recruitment and composition data on the associated herds. Collecting data from herds located through radio-tracking of collared cows minimized detection bias that might be present with opportunistically sampling any herd for composition and calf: cow ratios. Upon detection of collared cows a herd count was conducted and the animals were classified by sex and age (calves, yearlings, cows, or bulls). Sex was primarily determined using characteristics of the genital and abdominal region. Calves were not sexed. Yearling classification was not always possible near October given the similarity to subadults as the animals neared their third year. While smaller to medium-sized herds could be counted and classified accurately and completely (i.e., a total count), in most instances large herds were substantially more difficult, if not impossible, to age and sex in their entirety. However, these herds comprised the majority of the population, often differed substantially in composition, and could have different rates of calf survival from smaller herds (Estes, 2014). Thus, they were critical to estimates of calf recruitment and sex and age composition. To address this we developed a herd sampling method based on the position of the focal collared cow in the herd. Upon visual location of the collared cow in a large herd we classified a sample of animals in the section of the herd radiating out from around the collared cow until we were unable to accurately classify individuals by sex and age, at which point the observation stopped. Herds sampled in this way averaged 685 animals, of which an average of 57 animals were classified. While some biases by age and sex can be present in sampling herds, the loosely defined aggregations typical of large herds (not under immediate predation pressure) made other methods unfeasible given that herd size and composition changed continually and spacing precluded accurate and complete sexing and aging in large groups. Using collared cow locations as the sampling commencement point provided a degree of randomization by avoiding potential systematic biases in composition resulting from an animal's peripheral or interior herd position (Estes, 2014) and allowed us to accurately record data on herd composition and calf: cow ratios. A total of 6734 total counts and sample counts (hereafter herd counts) were conducted. We retained 5440 of these for fitting the countbased vital-rate model (see below), excluding counts near and during calving season (1 September to 31 October) when the presence or absence of a calf with a collared cow could not be reliably determined, and excluding counts during incomplete years at the start and end of the study. This yielded 241,251 records of whether an animal was a calf or a cow (and not a yearling or bull).
We determined whether a collared cow had a calf by a composite of interactive behaviors between the animals including nursing, grooming, and spacing (following). Vital rate estimates based on the calves of collared cows were derived from 216 cow-years observing 94 cows. We increased field efforts in three periods: (1) August-September, in order to facilitate estimation of first-year recruitment, (2) October-November, in order to facilitate estimation of fecundity, and (3) April-May, in order to facilitate estimation of wet-season calf survival.
We quantified causes of mortalities using three independent sources of information: (1) monitoring of collared cows, (2) opportunistic detection of mortalities within the population as a whole in the course of year-round work on wildebeest and carnivores, and (3) hunt follows as described in Droge et al. (2017a). Cow collars were equipped with mortality sensors, which became active after 12 h of no movement. When a mortality signal was detected, the collar was located and the area was assessed for signs of predators, poaching, disease, or starvation as per prior work on predation Droge et al., 2017a). Typically, given the high density of hyenas (see below), carcasses rarely had much of anything left besides the collar, making it difficult to definitively determine cause of death. Calf mortality was determined when the calf was no longer detected with the collared cow, except after July when calves and cows tend to begin dissociating as pregnant cows near parturition (Estes, 2014).
Through concurrent studies we were able to quantify the composition of the predator population, which provided more specificity about predator influences on the wildebeest population, including both predation and predation risk effects (Christianson et al., 2018;Creel, 2018;Creel et al., 2017Creel et al., , 2019Droge et al., 2017aDroge et al., , 2017bDroge et al., , 2019M'soka et al., 2016, 2017. Wildebeest are the primary prey for hyena, lion, and African wild dog, and an important prey for cheetah.
Due to depletion by poaching and conflict, the GLE carnivore guild is unusual because of its low numbers of lions and dominance of spot- were also present, with African wild dogs absent from the intensive study area from 2015 to 2019. Cheetah were stable to increasing in numbers during the study and ranged widely throughout the GLE.
Poaching was a suspected cause of mortality for a number of collared wildebeest. We removed these animals from the analyses of vital rates to avoid complications to do with varying levels of uncertainty in knowing whether or not they were actually poached and the possible biases that poaching would introduce into the estimation of environmentally determined influences on mortality. We examined the consequences of this policy and found that it had little effect on the estimation of basic vital rates and that it led to stronger evidence for the various effects on mortality.

| Model for survival and recruitment
We quantified wildebeest survival by fitting a single survival model to three different types of field measurements. The approach is summarized here, with mathematical details provided in the Appendix A.
The model was based around Siler's (1979) 5-parameter expression of mortality hazard (i.e., risk) as a constant modified by additional hazard components representing the immature and senescent phases of an organism's life. We incorporated additional modifications that modulated the hazard by season, location, and year. In the calf-survival analysis, we also represented calf dissociation from their mothers as a "hazard" in the mathematical sense.
The season effect was represented as a dichotomy between wet and dry seasons, with a parameter β S measuring the degree to which hazard is greater in the wet season. The location effect was represented as a single continuous value representing distance along the major axis of migration from an arbitrary point at the southeast limits of the migration (Figure 1), with a parameter β L measuring the degree to which hazard is greater with distance away from the southeastern corner of the study area. Year effects were represented by a vector β Y of parameters for each year indicating the degree to which that year was more or less hazardous than other years. Calf dissociation was represented by a parameter measuring the increase in dissociation "hazard" with time since 0.75 years of age.

| Bayesian parameter estimation and formal inference
We estimated model parameter distributions using Bayesian Monte Carlo (MC) methods based on Metropolis sampling (Givens & Hoeting, 2005;Metropolis et al., 1953;Smith, 2007). A specific likelihood function for the relevant model parameters was developed for each data set (collared cow survival, herd counts, and calf detection). Each model run involved six chains with at least 10,000 steps and adequate burn-in. We used a sequential Bayesian approach to unify the inference resulting from different model fits (Daniels et al., 2014;Garrard et al., 2012;Kurota et al., 2009;Michielsens et al., 2008). Different subsets of the model were fit to different data sets in sequence, with the posteriors from each fit becoming the priors for subsequent fits (see Appendix A).
Formal inference was derived directly from the estimated parameter distributions. These were interpreted using a variety of metrics-depending on the parameter in question-including expected values (EV), credible intervals (CI) based on highest posterior densities, one-sided lower credible limits (LCL0), and probabilities of being greater than zero (Pr > 0). Each of these metrics essentially describes a probability or quantile of probability that constitutes some form of evidence that can be compared to its complement in an "evidence ratio": ER = Pr/(1-Pr). We used the terms "substantial", "strong", "very strong", and "decisive" to interpret log10 evidence ratios (LERs) of 0.5, 1.0, 1.5, and 2.0, respectively (Kass & Raftery, 1995). Credible intervals were computed at the "strong" level of evidence, that is, 90.9% = 10÷(10 + 1).
We computed the distributions of derived parameters (e.g., recruitment, life expectancy, etc.) by applying their derivation to every element in the chain and computing metrics on the result.
To interpret interannual variation we computed both the range and the mean absolute difference (MAD) among all years of the interannual parameters, and we did this for every element in the chains. In the case of ratios among interannual hazard parameters, we used the geometric equivalents of the above metrics.
Specifically the maximum ratio (maxR) among any vector, h, of hazards was computed as exp(range(ln(h))), and the mean ratio (meanR) was computed as exp(MAD(ln(h))). For the purpose of calculation of this statistic, mean first-year hazard h1 for any given year was computed from first year survivorship S 1 as h 1 = − lnS 1 , and the interannual component of adult hazard was computed as hA = exp(β Y ).

| Longevity and life expectancy
We estimated longevity directly from modeled survivorship as the age attained by 0.1% of females. We estimated life expectancy at birth as the integral of the survivorship function (calculated numerically because there is no closed-form antiderivative of the Siler survivorship) (Canudas-Romo et al., 2018).

| Population trajectory
We examined the population trajectory both through direct aerial surveys of wildebeest, and by projection of vital rates using a demographic population model. The aerial survey data provided estimates of population size and hence an indication of whether there was evidence that the population was stable, increasing, or declining. The model projections explored the self-consistency of the estimated vital rates, and provided an alternate perspective on the question of population trend. If the model projections predicted a population explosion or crash with a severity that was inconsistent with the aerial survey data, this would be an indication of a problem with the vital-rate estimates or the assumptions behind them, or with the accuracy of the aerial survey estimates. Conversely, if the model projections were generally consistent with the aerial survey data, this would support the notion that the study had successfully documented several of the fundamental demographic characteristics of the population, including fecundity, recruitment, age-dependent survival, longevity, age structure, and sensitivity to season and spatial location.
Aerial population assessments were completed 12 times since 1970 using various methods. Partial-area surveys by incompletely documented methods were conducted in 1970, 1991, and 2001(Estes & East, 2009Kamweneshe et al., 2003;Tembo & Siawana, 1993). A more standardized and well-documented Systematic Reconnaissance Flight (SRF) sample methodology was used by Viljoen in 2004Viljoen in , 2007Viljoen in , 2009Viljoen in , 2011Viljoen in , 2013Viljoen in , and 2015 Jolly's II method to estimate population size and 95% confidence intervals (Viljoen, 2013). More recently, African Parks carried out total area counts (Norton-Griffiths, 1978, 2017 photogrammetry to improve the accuracy of counting large wildebeest herds (African Parks, 2018).
To explore the self-consistency of vital rates, we estimated their consequent population trajectories using a demographic population model based on standard Leslie matrix principles. A Leslie matrix model is comprised of a vector of abundances and a transition matrix. The abundance vector represented female wildebeest abundance in each of 51 six-month age classes, ranging from age 1 to 25 years (i.e., well above the plausible maximum).
The projection matrix was populated with adult survival rates from the adult survival model runs (model CC, see Table A1) and recruitment rates computed using fecundity and calf survival rates from the herd count model runs (model HC, see Table A1) and an assumed calf sex ratio of 0.5. Recruits were placed into the 1.0-year age class. Projection from one time-step to the next occurred over six-month intervals. Fecundity was assumed to be equal for all adults 2.0 years and older, and zero for yearlings. This is an approximation that maintains consistency with the implicit assumption of the herd composition analysis that only cows 2.0 years and older were reproductive. Absolute abundance values were irrelevant because the model did not incorporate density dependence. So, we initialized total abundance at an arbitrary value of 1 and referred to this as a "relative population size" in the results. Demographic population projection model simulation runs are typically initialized using either an observed age distribution or the stable age distribution (Lacy et al., 2021), the latter being the age distribution to which any standard Leslie matrix model will converge under constant vital rates, regardless of the initial distribution (Caughley, 1977;Sinclair et al., 2006). If both options were available, we would argue for use of the stable age distribution in the present situation because it represents a more neutral starting point for the exploration of the general self-consistency of vital rates estimated over an 8-year period-but the argument is moot because a sufficient observed age distribution was not available.
An observation age distribution could potentially be estimated from herd counts and tooth measurements during collaring, but the attempt would be hampered by under-representation of collared animals between the calf and full-grown adult stages and no individual collaring year exhibiting a large enough sample size across all age classes. Stable age distributions have been used by others in similar simulation contexts (Curtis et al., 2015;Monson et al., 2000;Taylor et al., 2009;Wiens, 2017). We computed the stable age distribution using a 50-year warm-up period at the start of each model run. During the warm-up period the total abundance was held at a constant value of 1 by dividing all values by their sum after each iteration. This was done in order to facilitate examination of convergence to stable age structure, and to avoid creating a misperception that absolute abundance was relevant to this population projection exercise. We confirmed that 50 years was sufficient time for convergence by verifying that there was effectively zero change in age distribution between the last two iterations. We could potentially have instead computed the stable age distribution as the dominant right eigenvector of an annualized transition matrix (Caughley, 1977;Sinclair et al., 2006), but since our model was biannual with vital rates dependent on season, this would have been arguably more complex than simply allowing the model to converge using a warm-up period as has been done by others in similar situations (Grear & Elderd, 2008;Wiens, 2017).
After the warm-up, the total abundance was allowed to fluctuate for five years and the annual population growth rate (λ) was computed as the one fifth power of the ratio between final and initial postwarm-up abundances. A five-year interval was chosen because we judged this to be the longest interval before unmodeled density dependent effects would need to be incorporated, and such effects were irrelevant to the scope of the exercise to examine self-consistency of vital rates during the study period.
The population projection was repeated for 5000 model runs with model parameters sampled from the posterior distributions of the adult survival and herd-count observational models in order to compute highest posterior density credible intervals for both λ and the stable age structure.

| RE SULTS
Results are summarized in Table A2  Wildebeest calving typically occurred in the northern-central portions of the park, with herds arriving on their wet season range several weeks after peak calving. The peak of the calving season was approximately October 1st.
Interannual variability was pronounced-the mean ratio in hazard between any two years was 1.53 (LCL 0 = 1.30, Model CC.Y) and the ratio in hazard between the highest and lowest years was 2.77

| Fecundity, calf survival, and recruitment based on herd-composition data
From herd composition data, mean fecundity was estimated to be 0.68 (0.57-0.82), mean survivorship to age 1.0 was estimated to be 0.56 (0.47-0.65), mean recruitment to age 1.0 was estimated to be 0.38 (0.36-0.40), and mean recruitment to age 2.0 was estimated to be 0.31 (0.27-0.36) (Model HC). There was strong evidence for pronounced interannual variation in each of these rates, that is, that mean absolute between-year differences in fecundity, survivorship, and recruitment were at least 0.11, 0.11, and 0.08, respectively (LCL 0 , Model HC.Y). Estimated rates for specific F I G U R E 2 Wildebeest survival metrics from birth to senescence estimated by a Siler hazard model fitted to herd-composition and collared-cow data (models HC and CC): (a) frequency (number of study periods) of ages of cows at mortality (dark lines) versus being observed alive (gray lines), (b) hazard, (c) instantaneous annualized survival rate computed directly from hazard, (d) survivorship. Dashes indicate 91% credible intervals. years are listed in Table A2. The degree of interannual variation was similar between calves and adults; the mean ratio in first-year calf hazard between any two calendar years was 1.64 (LCL 0 = 1.41, Model CC.Y), which is similar to the value of 1.53 reported for adults (above). There was no evidence for a decline in fecundity with cow age.

| Fecundity, calf survival, and recruitment based on calf-detection data
Mean fecundity before and after age 8 was 0.67 and 0.72, respectively (based on sample sizes of 189 and 60, respectively).

| Longevity and life expectancy
The estimated age exceeded by only 0.1% of females was 14.41 years , Figure 5a). The estimated life expectancy at birth was 3.71 years (CI 2.43-4.81).

F I G U R E 3
Wildebeest fecundity, survivorship, and calf dissociation from their mothers. Fecundity is plotted as the starting value for survivorship. Panel (a) compares estimates based on herd composition (HC) and collared-cow calf detections (CD). Panel (b) illustrates calf dissociation by plotting the probability that a calf is associated with a given cow on the same scale as estimated survivorship.

F I G U R E 4
Population size estimated from aerial surveys between 1970 and 2018. The single error bar in 1970 represents Berry's qualitative speculation as to how much the true population size may have been above their actual count. Error bars for the surveys carried out between 2009 and 2015 represent Viljoen's 95% confidence interval based on the Jolly's II method for sample counts. Surveys from 2016 to 2019 were total counts and hence did not quantify confidence intervals.

| Causes of death
We examined cause of death for 1086 carcasses (Figure 6). This included 981 animals with a known, natural cause of death, that is, where the cause of death could be reasonably ascertained and was not suspected of being unnatural (e.g., poaching). Of these, 67 were from collared cows, 215 were from hunt follows, and 699 were F I G U R E 5 Demographic population model predictions derived from estimated vital rates: (a) stable age structure after model warm-up, (b) population trajectory, with abundance represented as "relative population size" normalized to a value of 1 before each iteration of the warm-up period, and thereafter allowed to vary. The 50-year warm-up period achieved stabilization of the age distribution to provide a consistent starting point for estimation of population growth rate. F I G U R E 6 Estimated cause of death for 1086 documented wildebeest mortalities discovered by various means. "Birthing problems" included stillborn calves and mothers who died giving birth. "Calf starvation" includes calves that died after losing their mothers. Causes listed as "non-natural" include poaching, harvest for lion management purposes, entanglement in a fence, and encounter with a temporary electric fence used around a temporary airstrip in 2011-2012.  (Eberhardt, 2002); however, few aging techniques have been applied to African ungulates to allow precise age estimation of adults, and a novel age estimation technique from tooth allometry enabled evaluation of age-structured vital rates for wildebeest (Christianson et al., 2018). The estimated age exceeded by only 0.1% of females (14.4 years, CI 13.3-15.5) was somewhat consistent with the ages we obtained from tooth cementum annuli (oldest tooth aged: 11.1 years) and one to two years shorter than maximum longevity values reported from other populations: 17 years (Talbot & Talbot, 1963), 16 years (Watson, 1967, as cited by Estes, 2014), 16 years (Andere, 1975), 16-17 years (Attwell, 1982).
The estimated life expectancy at birth (3.7 years, CI 2.4 to 4.8) was shorter than other values in literature: 5.4 years (Andere, 1975), 7 years (Watson, 1967 as cited by Estes, 2014). Not surprisingly, wildebeest cows exhibited strong age-based survival rates, with older cows having lower rates of survival ( Figure 2). As with prior studies in the GLE, predation was strongly age-class specific, with hyenas and lions predating primarily adult wildebeest, and cheetah and wild dogs primarily calves and yearlings Droge et al., 2017a). The one exception to this was in the first few weeks of calving, when hyenas focused heavily on vulnerable newborn calves before resuming typical patterns of larger adult age classes; this was consistent with findings from the Serengeti-Mara wildebeest studies (Estes, 2014). Numerous studies have demonstrated age-based predation on adult ungulates, with coursing predators typically taking older animals in the adult age class in addition to the juvenile age class (Becker et al., 2009;Husseman et al., 2003;Kruuk, 1972;Kunkel et al., 1999;Schaller, 1972). Given spotted hyena (a coursing predator) were the dominant predator on wildebeest this would be expected, and the observed age dependency of mortality risk supported this (Figure 2b). While it was difficult to determine cause of death for calves of collared cows, concurrent lines of investigation indicated predation was again the main cause of mortality. Spatial data on wildebeest indicated that calving typically occurred in the northern portion of the LPNP, where predation risk was considerably lower, and calf: cow ratios quickly dropped in the first few weeks before slowing. Predation was most intense in the first few weeks following calving before hyena also resumed hunting primarily larger age classes (ZCP unpublished data). This supports identical findings from Estes (2014) indicating that hyena specialized on calves during the short window of the first few weeks of life, before predating larger age classes predominantly. Bolger et al. (2008) emphasized the need for evaluating factors limiting migratory ungulates across all phases of migration, given that these factors can vary with different phases, and our results support this recommendation. Adult wildebeest survival varied substantially between the wet and dry season range, primarily due to the significant increase in predation during the wet season. Herds spent the wet season in a relatively small area comprising the ecosystem's highest density of predators, before migrating north with the onset of the dry season into a much lower predator density.
Thus, migration to the more northern areas during the dry season and for calving is likely undertaken at least in part to escape the intense predation pressure of the wet season range. The determination of whether predators limit prey rests on several observations (Mills, 2007): the predation rate must be high, a substantial amount of the predation must be additive and not compensatory, and the killed prey must have high reproductive value. The evidence strongly supports each of these criteria. The predation rate is high; annual population mortality is high (estimated at 15.3%) and most of this is due to predation (81.4% or 98.0%, depending on how "possible" predations are accounted). If the predation was compensatory, we would expect evidence of starvation in the population, but starvation only accounted for 0.8% of known natural mortalities observed by unbiased means. Finally, high reproductive value individuals are predated, with the average age at mortality of predated collared females being 7.80 years and no evidence for a decline in fecundity this could be as the result of wildebeest selecting areas of poorer habitat quality to avoid predictable areas of high predation risk, consistent with the Control of Risk Hypothesis (Creel, 2018).
While predation and predation risk effects exerted a strong effect on wildebeest demography and likely migration, we cannot rule out the importance of other drivers in seasonal movements. Ecological factors such as access to forage and its component variables such as precipitation, fire, surface water availability, and nutrients have been identified as key drivers for migration for wildebeest and other migratory herbivores (Geremia et al., 2014;Holdo et al., 2009aHoldo et al., , 2009bHopcraft et al., 2015;Sinclair & Norton Griffiths, 1982). It is possible that the indigestible content of wet season forage (including water content) was so high that animals could not easily maintain nutrition (Hopcraft et al., 2001(Hopcraft et al., , 2010. Given the extensive flooding of the GLE in the rains-and subsequent high water-content of forage-it is possible that nutrient limitation drove wildebeest herds to occupy areas with high predation risk in order to fulfill energy, protein, and possibly mineral requirements. While it is not the focus of this paper, future analyses will be focused on the drivers of wildebeest spatial dynamics and migration in the GLE. Poaching of wildebeest also impacted wildebeest survival in the ulation to decline more than 20%. This is by no means assured, given the disparate investigations underlying these estimates-including tooth allometry, herd counts, and traditional mark-resight survival estimation-and the absence of any stabilizing density-dependence mechanism built into the model. It also supports the notion apparent from the aerial population survey data that there is no clear evidence of upward or downward trend in GLE population size at present, notwithstanding the considerable imprecision brought about by large sample errors in SRF sampling and potential for observer bias through over-or underestimation (particularly of very large herds of animals) inherent in aerial surveys prior to the use of photogrammetry (Lamprey et al., 2020;Schuette et al., 2018). A more-detailed future analysis could consider both non-zero yearling fecundity, and 2-year-old fecundity being lower than that of older adults-both of which have been noted in the literature with substantial variability between years and between systems (Mduma et al., 1999;Mills & Shenk, 1992;Tambling & Du Toit, 2005), and for which we have a limited number of supporting observations from the GLE.
Nevertheless, the increased resource protection efforts for both the GLE and its wildlife since 2003 has resulted in very high rates of survival for wildebeest-dependent carnivores such as hyena (M'soka et al., 2016), and no evidence of population declines, indicating carnivore populations are likely tracking herbivore populations that are stable to increasing. The influence of predation on wildebeest demography relative to resource-limitation was in striking contrast to studies elsewhere such as that of Mduma et al. (1999) which found forage to be the primary limiting factor for Serengeti wildebeest.
Whether the strong predation limitation is a factor of current wildebeest population size or transitional recovery dynamics, or whether it reflects differences in ecological dynamics in the GLE versus other populations remains to be seen. Similarly, dynamics of wildebeest could change significantly with changes in diet and preference of hyenas and other large carnivores. Increases in predation on other relatively abundant species such as zebra could either exert a stabilizing influence on wildebeest through prey switching or increase predation impacts through apparent competition. Continued restoration of both the migrant and resident herbivore community in the GLE will likely have an array of impacts on predator-prey dynamics .
Large herbivore migrations are in decline worldwide and while Africa has the largest number of remaining migrations, only one-the Serengeti-Mara Ecosystem-is mostly protected (Harris et al., 2009), although even this range is at risk from road developments and other human impacts (Holdo et al., 2011;Hopcraft et al., 2015). Harris et al. (2009) identified key actions needed to conserve migratory populations, specifically the need to secure seasonal range, resource protection, governmental support, and minimizing fences. Migratory populations are typically more abundant than resident ones as they can likely avoid the carrying-capacity restraints imposed on residents by tracking resources and avoiding predation; however, when constraints or barriers to migration occur, dramatic declines and even collapses in populations can occur (Bekenov et al., 1998;Berger, 2004;Berry, 1997;East, 1988;Estes & East, 2009;Whyte & Joubert, 1988), and wildebeest migrations appear particularly sensitive to disruption (Holdo et al., 2011;Morrison et al., 2016;Morrison & Bolger, 2014;Thirgood et al., 2004).
The GLE wildebeest population is strongly limited by predation, particularly on their wet season range that encompasses the system's highest densities of predators. Rapid human encroachment in the form of agricultural conversion in the UWZGMA that comprises much of the GLE wildebeest dry season range and lies completely outside the LPNP is threatening both the current migratory range of the wildebeest and any potential expansion and resumption of historical transfrontier migrations. If wildebeest were to be decoupled from their migratory range outside the LPNP then the population would likely be subjected to intense year-round predation from a growing carnivore population in the south, which would likely drive declines if herds were not able to escape high levels of predation on their dry season range and during calving. Consequently, there is an urgent need for protection of the entire GLE wildebeest range in order to secure the long-term viability of this keystone species in a rapidly changing ecosystem. Increased protection of current and potential dry season range outside LPNP should be considered the highest priority for conservation of wildebeest and other dependent species .

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interests in relation to the publication of this paper.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data supporting the results of this paper are archived at Dryad at https://doi.org/10.5061/dryad.0k6djhb3f.

AG E O F CO LL A R E D COWS
We estimated the age of collared cows using tooth measurements in a two-stage process based on the cementum annuli of the first and second incisors (Christianson et al., 2018). The first stage utilized teeth extracted from adult wildebeest carcasses detected in the course of fieldwork on predators and wildebeest. These teeth were aged at the Matson Lab (Manhattan, MT, USA) using cementum annuli techniques, from which we developed relationships between these age estimates and the physical dimensions of the teeth. Teeth from found carcasses ranged in age from 1.1 to 11.1 years, based directly on cementum annuli measurement. The second stage utilized tooth measurements taken from immobilized wildebeest during collaring and applied the carcass-based relationships to yield an estimated age for each collared wildebeest. We made slight adjustments to the lab-reported ages to accommodate the constraint that all wildebeest in the study system are born at a certain time of year (October). We measured three tooth dimensions: arcade breadth (AB), labio-lingual width (LLW), and emergent crown height (ECH) (Christianson et al., 2018). Each of these can be measured both with a carcass tooth in hand, and from the in-situ teeth of an immobilized animal. We regressed cementum age against each of the three dimensions in turn ( Figure A1) and found that AB explained the most variation (R 2 = 0.68, N = 61, p < .0001), followed by ECH (R 2 = 0.59, N = 61, p < .0001) then LLW (R 2 = 0.41, N = 61, p < .0001). We, therefore, used AB as the primary basis for age estimation for collared animals, except where there were data gaps in AB, in which case we used ECH, and then LLW, as necessary. A reviewer suggested combining all three metrics in order to achieve a more accurate estimate.
This would be worth exploring more, although we note that certain multimetric approaches would be precluded by data gaps. Some animals were recollared after a few years, in which case we used the age estimate from the first collaring because these were generally more precise. Age estimates from recollaring generally agreed with those from initial collaring but in some cases were relative underestimates.

F I G U R E A 1
Estimation of wildebeest age (years) from the dimensions of teeth obtained from found carcasses. The dimensions for collared wildebeest at the time of collaring are also given, in order to show that age estimation for collared animals did not rely on excessive extrapolation beyond the range of measurements available from carcasses.
Finally, we qualitatively compared the distribution of tooth dimensions between carcasses and collared animals to verify that the carcass-derived teeth were representative and that age-estimation for collared animals would not rely on excessive extrapolation.

M O D E L FO R S U RV I VA L A N D R ECRU ITM E NT
We quantified wildebeest survival by fitting a single survival model to three different types of field measurements, as described below.
The model was of the form introduced by Siler (1979), whereby total hazard, h, is the summation of terms for (1) immaturity hazard, (2) constant hazard, and (3) senescence hazard: where t represents age in years, and a 1 , b 1 , a 2 , a 3 , b 3 are parameters with subscripts representing the three forms of hazards listed above.
We use the terms "hazard" and "mortality risk" interchangeably depending on context, the former being more prominent in mathematical survival analysis.
The hazard function leads to a function for survivorship, S (Lee & Wang, 2003): In the original Siler model, parameters a 1 , a 2 , and a 3 specify the magnitude of each hazard, and parameters b 1 and b 3 specify the rate of reduction in immaturity hazard and increase in senescence hazard, respectively. We expressed a 1 and a 3 as a function of more biologically meaningful parameters S 1 and t f , where S 1 specified the first-year survivorship, and t f specified the time until survivorship fell to a specified low level, f = 0.01 (Barlow & Boveng, 1991;Becker et al., 2013). Given S 1 and t f , as well as b 1 , a 2 , and b 3 , values of a 1 and a 3 can be calculated according to re-arrangement of the original survivorship function: These two equations are interdependent, but their solution converges in two iterations under typical values of the known parameters.
We explored models with several subsets and variants of the survivorship function ( To explore spatial variation in hazard, a location effect, X L , was represented with a fixed parameter X L,ref denoting a reference location, and a fitted parameter β L : Location was represented as the relative distance, X L ∈ [0,1], along the major axis of wildebeest migration measured from an arbitrary point at the southeast limits of the migration. The reference location was chosen as the median relative distance, X Lref = 0.355.
All secondary effects were centered around zero and then exponentiated to prevent negative hazards as indicated in the above equations.
A calf dissociation effect was represented as a "hazard." Calves begin to dissociate from their mothers at approximately 0.75 years of age (Estes, 2014). This has the same numerical effect as mortality

M O D E L FIT TI N G TO K N OWN FATE S O F A D U LT S
We estimated adult survival using a known-fates approach, assuming perfect detection of living animals and their deaths. We assured perfect detection through frequent attempts to relocate each collared animal, rotating through all animals and focusing on the ones seen least recently, and grouping the resulting relocation data into six-month periods. Six months is longer than the time it took to determine the state of all collared animals, given the level of year-round search effort we maintained (800-1100 person-days per year). We assumed a closed population. No other populations are nearby and opportunity for emigration is strongly limited by landscape features.
Four collars were lost due to either malfunction or destruction by predators, scavengers, or poachers. In these cases, we removed the animal from the study at the start of the period in which the collar loss occurred. Thus, the animal had a known fate (survival) through to the period prior to the collar loss. We encountered no evidence that collar loss was correlated with fate or predictor variables such as season, location, or age.
The likelihood function for observation of survivorship over a known period of time to a known fate (living, or died at a known time) is Rodriguez (2010): where t i,j is the amount of time (years) in which wildebeest I was in the study during semester j. This value is 0.5 years for wildebeest who were collared before a semester and survived past the end of that semester, but <0.5 years for wildebeest either collared or determined to have died during a semester (or both), and 0 in semesters before collaring or after death.

M O D E L FIT TI N G TO O B S E RVATI O N S O F H E R D CO M P OS ITI O N
We estimated recruitment from fecundity and calf survival in two ways: (1) from calf: cow ratios in herd counts, and (2) from intensive monitoring of the calves of collared cows, thereby providing two concurrent and independent measures . Animals classified as calves or cows in herd counts can be viewed as an outcome of a binomial process with a probability (θ) that each calf or cow (Y) in the herd count will be either a calf (Y = 1) or a cow (Y = 0): A standard binomial likelihood accompanies such a model: where Θ is a vector of all model parameters, y is the number of calves, and n is the number of calves and cows. These likelihoods were raised to a power of a 0.2 to compensate for the estimated rate of repeat records of the same animal in a given year. Recruitment was not a directly fitted model parameter. Rather, it was a model output computed from other modeled quantities depending on the fitting context-either from F and S under the calf detection data, or from θ under herd composition data (see below). Recruitment to specific ages 1 and 2 was denoted R 1 and R 2 , respectively.

M O D E L FIT TI N G TO D E TEC TI O N S O F TH E C A LV E S O F CO LL A R E D COWS
The second way in which we estimated recruitment from fecundity and calf survival was to monitor the association of calves with individual cows, that is, by recording at each resighting of a collared cow the detection or nondetection of a calf associated with that cow.
Calves stay with their mother for approximately 9 months, at which time they can begin to dissociate from their mother, until all are dissociated once new calves arrive. Twins and adoption do not occur (Estes, 2014, ZCP unpublished data) and we assumed that calves died if their mother died. With a few minutes of observation, it is generally straightforward to determine if a cow has a calf and to which cow a given calf belongs; we referred to this as a "true detection".
Such determinations are not infallible, however, and it possible to erroneously record another cow's calf as being associated with a specific collared cow. We referred to this as a "false detection". We accounted for this with an observational model that incorporated true and false detection probabilities.
The collared-cow calf detection model is not unlike traditional capture-resight models, but the possibility of false detections necessitated special treatment via a likelihood function that accounted for the specific observational circumstances of our study using standard binomial and multinomial likelihood principles: where F is fecundity, q is the false detection probability, Y is a ragged where z indexes the last occasion when the calf was still alive. For example, the detectable life history for a calf that was born but never detected (z = 0) would be {0,0,…,0}; the history for a calf that died after the third sampling occasion (z = 3) would be {1,1,1,0,0,…,0}; and the Z n i ,z = 1 1 , 1 2 , ⋯ , 1 z , 0 1 , 0 2 , ⋯ , 0 n i −z , for z ∈ 0, ⋯ , n i , with the shorthand Z = Z n i ,z history for a calf that survived at least through all sampling occasions (z = n i ) would be {1,1,…,1}. A multinomial likelihood is summed over all of these possibilities, weighting each one by the appropriate interval of survivorship, (ΔS) i,z , between consecutive sampling occasions (Choquet et al., 2011;Fieberg & DelGuidice, 2011;Friedman, 1982): which leads to an overall fecund likelihood: where p is the true detection probability. The six individual summed terms in the fecund likelihood explore the six possible ways in which each detection (Y i,j ∈ {1, 0}), can arise from the calf being alive or not (Z j ∈ {1, 0}), the correct calf being detected (p), or another calf being detected (q) (given that the cow was fecund in that calf-year). The barren fecundity is simpler, needing only to account for two detection possibilities, and no calf life history: At the outset, we considered that the two methods of estimating recruitment should lead to similar estimates, but they each may have different strengths. We considered herd counts to provide the best overall estimates of recruitment, being based on a very large sample size. Monitoring of the calves of individual cows is based on a much smaller sample size but is the only approach that allows incorporation of covariates that rely on movement history. It also provides a useful independent check on the herd-count approach.

M O D E L FIT TI N G S EQ U E N CE
The model was fit to the different data sets using a sequential Bayesian approach, with the posterior parameter distributions from initial fits informing the prior parameter distributions of subsequent model fits (Daniels et al., 2014;Garrard et al., 2012;Kurota et al., 2009;Michielsens et al., 2008). The sequence was initiated by fitting the adult survival model to known-fates data with noninformative priors and no inter-annual variation (Model CC0, Table A1, Figure A2). This allowed some specification of the constant hazard (a2) parameter, which was found to be a necessary prior for convergence of the herd-counts models (Models HC and HC.Y). Fitting of the HC model provided some specification of the first-year survivorship (S1) which enable more precise fits of the known-fates models (Models CC, CC.S, CC.L, and CC.Y) and estimation of associated adult survival parameters (βS, βL, βY). Finally, the calf-detection models (CD and CD.Y) used posteriors from the HC.Y model as priors for the rate of reduction in immaturity hazard (b1) and constant hazard (a2) in order to achieve convergence, leaving the primary recruitment parameters (R, F, and S1) free to be estimated from noninformative priors. As an alternative to this sequential Bayesian approach, a single fit to an integral model with all data sets might be possible, but as Staton et al. (2017) demonstrated with a similar example in fisheries population dynamics, the results would likely be similar and the outward simplicity of an integrated approach might not outweigh the underlying computational expense. The sequential approach is computationally quicker and more assured of convergence because it explores a sequence of small low-dimensional parameter spaces instead of one vast high-dimensional parameter space. Our final sequence took about 7 h to run on a personal computer with a 2.8 GHz 4-core processor, and model development involved hundreds of test runs.   Table A1 for explanation of terms. Parameter estimation relied on a sequential Bayesian approach in order to estimate distributions of parameters that are determinant of both calf survival and adult survival and thus informed by multiple types of observation.